Stokes I Simultaneous Image and Instrument Modeling

In this tutorial, we will create a preliminary reconstruction of the 2017 M87 data on April 6 by simultaneously creating an image and model for the instrument. By instrument model, we mean something akin to self-calibration in traditional VLBI imaging terminology. However, unlike traditional self-cal, we will at each point in our parameter space effectively explore the possible self-cal solutions. This will allow us to constrain and marginalize over the instrument effects, such as time variable gains.

To get started we load Comrade.

using Comrade
  Activating project at `~/work/Comrade.jl/Comrade.jl/examples`
using Pyehtim
using LinearAlgebra

For reproducibility we use a stable random number genreator

using StableRNGs
rng = StableRNG(123)
StableRNGs.LehmerRNG(state=0x000000000000000000000000000000f7)

Load the Data

To download the data visit https://doi.org/10.25739/g85n-f134 First we will load our data:

obs = ehtim.obsdata.load_uvfits(joinpath(dirname(pathof(Comrade)), "..", "examples", "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))
Python: <ehtim.obsdata.Obsdata object at 0x7fe8b1d6d9c0>

Now we do some minor preprocessing:

  • Scan average the data since the data have been preprocessed so that the gain phases coherent.
  • Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.
obs = scan_average(obs.add_fractional_noise(0.02))
Python: <ehtim.obsdata.Obsdata object at 0x7fe8b1d43f10>

Now we extract our complex visibilities.

dvis = extract_table(obs, ComplexVisibilities())
EHTObservation{Float64,Comrade.EHTVisibilityDatum{Float64}, ...}
  source: M87
  mjd: 57849
  frequency: 2.29070703125e11
  bandwidth: 1.856e9
  stations: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples: 284

##Building the Model/Posterior

Now, we must build our intensity/visibility model. That is, the model that takes in a named tuple of parameters and perhaps some metadata required to construct the model. For our model, we will use a raster or ContinuousImage for our image model. The model is given below:

function sky(θ, metadata)
    (;fg, c, σimg) = θ
    (;ftot, K, meanpr, grid, cache) = metadata
    # Transform to the log-ratio pixel fluxes
    cp = meanpr .+ σimg.*c.params
    # Transform to image space
    rast = (ftot*(1-fg))*K(to_simplex(CenteredLR(), cp))
    img = IntensityMap(rast, grid)
    m = ContinuousImage(img, cache)
    # Add a large-scale gaussian to deal with the over-resolved mas flux
    g = modify(Gaussian(), Stretch(μas2rad(250.0), μas2rad(250.0)), Renormalize(ftot*fg))
    return m + g
end
sky (generic function with 1 method)

Unlike other imaging examples (e.g., Imaging a Black Hole using only Closure Quantities) we also need to include a model for the instrument, i.e., gains as well. The gains will be broken into two components

  • Gain amplitudes which are typically known to 10-20%, except for LMT, which has amplitudes closer to 50-100%.
  • Gain phases which are more difficult to constrain and can shift rapidly.
function instrument(θ, metadata)
    (; lgamp, gphase) = θ
    (; gcache, gcachep) = metadata
    # Now form our instrument model
    gvis = exp.(lgamp)
    gphase = exp.(1im.*gphase)
    jgamp = jonesStokes(gvis, gcache)
    jgphase = jonesStokes(gphase, gcachep)
    return JonesModel(jgamp*jgphase)
end
instrument (generic function with 1 method)

The model construction is very similar to Imaging a Black Hole using only Closure Quantities, except we include a large scale gaussian since we want to model the zero baselines. For more information about the image model please read the closure-only example. Let's discuss the instrument model Comrade.JonesModel. Thanks to the EHT pre-calibration, the gains are stable over scans. Therefore, we can model the gains on a scan-by-scan basis. To form the instrument model, we need our

  1. Our (log) gain amplitudes and phases are given below by lgamp and gphase
  2. Our function or cache that maps the gains from a list to the stations they impact gcache.
  3. The set of Comrade.JonesPairs produced by jonesStokes

These three ingredients then specify our instrument model. The instrument model can then be combined with our image model cimg to form the total JonesModel.

Now, let's set up our image model. The EHT's nominal resolution is 20-25 μas. Additionally, the EHT is not very sensitive to a larger field of view. Typically 60-80 μas is enough to describe the compact flux of M87. Given this, we only need to use a small number of pixels to describe our image.

npix = 32
fovx = μas2rad(150.0)
fovy = μas2rad(150.0)
7.27220521664304e-10

Now let's form our cache's. First, we have our usual image cache which is needed to numerically compute the visibilities.

grid = imagepixels(fovx, fovy, npix, npix)
buffer = IntensityMap(zeros(npix, npix), grid)
cache = create_cache(NFFTAlg(dvis), buffer, BSplinePulse{3}())
VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im  …  0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

Second, we now construct our instrument model cache. This tells us how to map from the gains to the model visibilities. However, to construct this map, we also need to specify the observation segmentation over which we expect the gains to change. This is specified in the second argument to jonescache, and currently, there are two options

  • FixedSeg(val): Fixes the corruption to the value val for all time. This is usefule for reference stations
  • ScanSeg(): which forces the corruptions to only change from scan-to-scan
  • TrackSeg(): which forces the corruptions to be constant over a night's observation

For this work, we use the scan segmentation for the gain amplitudes since that is roughly the timescale we expect them to vary. For the phases we use a station specific scheme where we set AA to be fixed to unit gain because it will function as a reference station.

gcache = jonescache(dvis, ScanSeg())
gcachep = jonescache(dvis, ScanSeg(); autoref=SEFDReference((complex(1.0))))

using VLBIImagePriors

Now we need to specify our image prior. For this work we will use a Gaussian Markov Random field prior Since we are using a Gaussian Markov random field prior we need to first specify our mean image. This behaves somewhat similary to a entropy regularizer in that it will start with an initial guess for the image structure. For this tutorial we will use a a symmetric Gaussian with a FWHM of 60 μas

fwhmfac = 2*sqrt(2*log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(50.0)./fwhmfac))
imgpr = intensitymap(mpr, grid)
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 NamedDimsArray(::Matrix{Float64}, (:X, :Y)):
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8
 (-3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
 (-3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
 (-2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
 (-2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
 (-2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46586e-6        5.12225e-6     2.46586e-6
  (2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
  (2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
  (2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
  (3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
  (3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
  (3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8

Now since we are actually modeling our image on the simplex we need to ensure that our mean image has unit flux

imgpr ./= flux(imgpr)
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 Matrix{Float64}:
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.38049e-8        1.3254e-7      6.38049e-8
 (-3.29522e-10)     1.3254e-7         2.75321e-7     1.3254e-7
 (-3.06796e-10)     2.62224e-7        5.4471e-7      2.62224e-7
 (-2.84071e-10)     4.94121e-7        1.02642e-6     4.94121e-7
 (-2.61345e-10)     8.86803e-7   …    1.84213e-6     8.86803e-7
 (-2.38619e-10)     1.51585e-6        3.14882e-6     1.51585e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46784e-6        5.12636e-6     2.46784e-6
  (2.38619e-10)     1.51585e-6        3.14882e-6     1.51585e-6
  (2.61345e-10)     8.86803e-7   …    1.84213e-6     8.86803e-7
  (2.84071e-10)     4.94121e-7        1.02642e-6     4.94121e-7
  (3.06796e-10)     2.62224e-7        5.4471e-7      2.62224e-7
  (3.29522e-10)     1.3254e-7         2.75321e-7     1.3254e-7
  (3.52247e-10)     6.38049e-8        1.3254e-7      6.38049e-8

and since our prior is not on the simplex we need to convert it to unconstrained or real space.

meanpr = to_real(CenteredLR(), Comrade.baseimage(imgpr))
32×32 Matrix{Float64}:
 -7.55422  -6.82317  -6.14085  -5.50727   …  -6.14085  -6.82317  -7.55422
 -6.82317  -6.09211  -5.4098   -4.77622      -5.4098   -6.09211  -6.82317
 -6.14085  -5.4098   -4.72748  -4.0939       -4.72748  -5.4098   -6.14085
 -5.50727  -4.77622  -4.0939   -3.46032      -4.0939   -4.77622  -5.50727
 -4.92243  -4.19137  -3.50906  -2.87548      -3.50906  -4.19137  -4.92243
 -4.38632  -3.65527  -2.97295  -2.33937   …  -2.97295  -3.65527  -4.38632
 -3.89895  -3.1679   -2.48558  -1.852        -2.48558  -3.1679   -3.89895
 -3.46032  -2.72927  -2.04695  -1.41337      -2.04695  -2.72927  -3.46032
 -3.07043  -2.33937  -1.65705  -1.02348      -1.65705  -2.33937  -3.07043
 -2.72927  -1.99821  -1.3159   -0.682317     -1.3159   -1.99821  -2.72927
  ⋮                                       ⋱             ⋮        
 -3.07043  -2.33937  -1.65705  -1.02348      -1.65705  -2.33937  -3.07043
 -3.46032  -2.72927  -2.04695  -1.41337      -2.04695  -2.72927  -3.46032
 -3.89895  -3.1679   -2.48558  -1.852     …  -2.48558  -3.1679   -3.89895
 -4.38632  -3.65527  -2.97295  -2.33937      -2.97295  -3.65527  -4.38632
 -4.92243  -4.19137  -3.50906  -2.87548      -3.50906  -4.19137  -4.92243
 -5.50727  -4.77622  -4.0939   -3.46032      -4.0939   -4.77622  -5.50727
 -6.14085  -5.4098   -4.72748  -4.0939       -4.72748  -5.4098   -6.14085
 -6.82317  -6.09211  -5.4098   -4.77622   …  -5.4098   -6.09211  -6.82317
 -7.55422  -6.82317  -6.14085  -5.50727      -6.14085  -6.82317  -7.55422

Now we can form our metadata we need to fully define our model.

metadata = (;ftot=1.1, K=CenterImage(imgpr), meanpr, grid, cache, gcache, gcachep)
(ftot = 1.1, K = VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}([0.9944957386363635 -0.0053267045454545746 … 0.005326704545454565 0.005504261363636378; -0.0053267045454545746 0.9948393969941348 … 0.005160603005865086 0.005326704545454551; … ; 0.005326704545454565 0.005160603005865086 … 0.9948393969941348 -0.005326704545454546; 0.005504261363636378 0.005326704545454551 … -0.005326704545454546 0.9944957386363636], (32, 32)), meanpr = [-7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; … ; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; -7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378], grid = GriddedKeys{(:X, :Y)}
	X: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
	Y: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
, cache = VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im  …  0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]), gcache = JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25)), gcachep = JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))]))

We will also fix the total flux to be the observed value 1.1. This is because total flux is degenerate with a global shift in the gain amplitudes making the problem degenerate. To fix this we use the observed total flux as our value.

Moving onto our prior, we first focus on the instrument model priors. Each station requires its own prior on both the amplitudes and phases. For the amplitudes we assume that the gains are apriori well calibrated around unit gains (or 0 log gain amplitudes) which corresponds to no instrument corruption. The gain dispersion is then set to 10% for all stations except LMT, representing that we expect 10% deviations from scan-to-scan. For LMT we let the prior expand to 100% due to the known pointing issues LMT had in 2017.

using Distributions
using DistributionsAD
distamp = station_tuple(dvis, Normal(0.0, 0.1); LM = Normal(1.0))
(AA = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AP = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AZ = Distributions.Normal{Float64}(μ=0.0, σ=0.1), JC = Distributions.Normal{Float64}(μ=0.0, σ=0.1), LM = Distributions.Normal{Float64}(μ=1.0, σ=1.0), PV = Distributions.Normal{Float64}(μ=0.0, σ=0.1), SM = Distributions.Normal{Float64}(μ=0.0, σ=0.1))

For the phases, as mentioned above, we will use a segmented gain prior. This means that rather than the parameters being directly the gains, we fit the first gain for each site, and then the other parameters are the segmented gains compared to the previous time. To model this we break the gain phase prior into two parts. The first is the prior for the first observing timestamp of each site, distphase0, and the second is the prior for segmented gain ϵₜ from time i to i+1, given by distphase. For the EHT, we are dealing with pre-2*rand(rng, ndim) .- 1.5calibrated data, so often, the gain phase jumps from scan to scan are minor. As such, we can put a more informative prior on distphase.

Warning

We use AA (ALMA) as a reference station so we do not have to specify a gain prior for it.

distphase = station_tuple(dvis, DiagonalVonMises(0.0, inv(π^2)))
(AA = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AP = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AZ = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), JC = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), LM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), PV = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), SM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688))

In addition we want a reasonable guess for what the resolution of our image should be. For radio astronomy this is given by roughly the longest baseline in the image. To put this into pixel space we then divide by the pixel size.

beam = beamsize(dvis)
rat = (beam/(step(grid.X)))
5.279832635689346

To make the Gaussian Markov random field efficient we first precompute a bunch of quantities that allow us to scale things linearly with the number of image pixels. This drastically improves the usual N^3 scaling you get from usual Gaussian Processes.

crcache = MarkovRandomFieldCache(meanpr)
VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}(sparse([1, 2, 32, 33, 993, 1, 2, 3, 34, 994  …  31, 991, 1022, 1023, 1024, 32, 992, 993, 1023, 1024], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 65, 65, 65, 65, 65, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 84, 84, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86, 86, 86, 86, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 101, 101, 101, 101, 101, 102, 102, 102, 102, 102, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104, 105, 105, 105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110, 111, 111, 111, 111, 111, 112, 112, 112, 112, 112, 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 116, 116, 116, 116, 116, 117, 117, 117, 117, 117, 118, 118, 118, 118, 118, 119, 119, 119, 119, 119, 120, 120, 120, 120, 120, 121, 121, 121, 121, 121, 122, 122, 122, 122, 122, 123, 123, 123, 123, 123, 124, 124, 124, 124, 124, 125, 125, 125, 125, 125, 126, 126, 126, 126, 126, 127, 127, 127, 127, 127, 128, 128, 128, 128, 128, 129, 129, 129, 129, 129, 130, 130, 130, 130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 132, 132, 133, 133, 133, 133, 133, 134, 134, 134, 134, 134, 135, 135, 135, 135, 135, 136, 136, 136, 136, 136, 137, 137, 137, 137, 137, 138, 138, 138, 138, 138, 139, 139, 139, 139, 139, 140, 140, 140, 140, 140, 141, 141, 141, 141, 141, 142, 142, 142, 142, 142, 143, 143, 143, 143, 143, 144, 144, 144, 144, 144, 145, 145, 145, 145, 145, 146, 146, 146, 146, 146, 147, 147, 147, 147, 147, 148, 148, 148, 148, 148, 149, 149, 149, 149, 149, 150, 150, 150, 150, 150, 151, 151, 151, 151, 151, 152, 152, 152, 152, 152, 153, 153, 153, 153, 153, 154, 154, 154, 154, 154, 155, 155, 155, 155, 155, 156, 156, 156, 156, 156, 157, 157, 157, 157, 157, 158, 158, 158, 158, 158, 159, 159, 159, 159, 159, 160, 160, 160, 160, 160, 161, 161, 161, 161, 161, 162, 162, 162, 162, 162, 163, 163, 163, 163, 163, 164, 164, 164, 164, 164, 165, 165, 165, 165, 165, 166, 166, 166, 166, 166, 167, 167, 167, 167, 167, 168, 168, 168, 168, 168, 169, 169, 169, 169, 169, 170, 170, 170, 170, 170, 171, 171, 171, 171, 171, 172, 172, 172, 172, 172, 173, 173, 173, 173, 173, 174, 174, 174, 174, 174, 175, 175, 175, 175, 175, 176, 176, 176, 176, 176, 177, 177, 177, 177, 177, 178, 178, 178, 178, 178, 179, 179, 179, 179, 179, 180, 180, 180, 180, 180, 181, 181, 181, 181, 181, 182, 182, 182, 182, 182, 183, 183, 183, 183, 183, 184, 184, 184, 184, 184, 185, 185, 185, 185, 185, 186, 186, 186, 186, 186, 187, 187, 187, 187, 187, 188, 188, 188, 188, 188, 189, 189, 189, 189, 189, 190, 190, 190, 190, 190, 191, 191, 191, 191, 191, 192, 192, 192, 192, 192, 193, 193, 193, 193, 193, 194, 194, 194, 194, 194, 195, 195, 195, 195, 195, 196, 196, 196, 196, 196, 197, 197, 197, 197, 197, 198, 198, 198, 198, 198, 199, 199, 199, 199, 199, 200, 200, 200, 200, 200, 201, 201, 201, 201, 201, 202, 202, 202, 202, 202, 203, 203, 203, 203, 203, 204, 204, 204, 204, 204, 205, 205, 205, 205, 205, 206, 206, 206, 206, 206, 207, 207, 207, 207, 207, 208, 208, 208, 208, 208, 209, 209, 209, 209, 209, 210, 210, 210, 210, 210, 211, 211, 211, 211, 211, 212, 212, 212, 212, 212, 213, 213, 213, 213, 213, 214, 214, 214, 214, 214, 215, 215, 215, 215, 215, 216, 216, 216, 216, 216, 217, 217, 217, 217, 217, 218, 218, 218, 218, 218, 219, 219, 219, 219, 219, 220, 220, 220, 220, 220, 221, 221, 221, 221, 221, 222, 222, 222, 222, 222, 223, 223, 223, 223, 223, 224, 224, 224, 224, 224, 225, 225, 225, 225, 225, 226, 226, 226, 226, 226, 227, 227, 227, 227, 227, 228, 228, 228, 228, 228, 229, 229, 229, 229, 229, 230, 230, 230, 230, 230, 231, 231, 231, 231, 231, 232, 232, 232, 232, 232, 233, 233, 233, 233, 233, 234, 234, 234, 234, 234, 235, 235, 235, 235, 235, 236, 236, 236, 236, 236, 237, 237, 237, 237, 237, 238, 238, 238, 238, 238, 239, 239, 239, 239, 239, 240, 240, 240, 240, 240, 241, 241, 241, 241, 241, 242, 242, 242, 242, 242, 243, 243, 243, 243, 243, 244, 244, 244, 244, 244, 245, 245, 245, 245, 245, 246, 246, 246, 246, 246, 247, 247, 247, 247, 247, 248, 248, 248, 248, 248, 249, 249, 249, 249, 249, 250, 250, 250, 250, 250, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252, 253, 253, 253, 253, 253, 254, 254, 254, 254, 254, 255, 255, 255, 255, 255, 256, 256, 256, 256, 256, 257, 257, 257, 257, 257, 258, 258, 258, 258, 258, 259, 259, 259, 259, 259, 260, 260, 260, 260, 260, 261, 261, 261, 261, 261, 262, 262, 262, 262, 262, 263, 263, 263, 263, 263, 264, 264, 264, 264, 264, 265, 265, 265, 265, 265, 266, 266, 266, 266, 266, 267, 267, 267, 267, 267, 268, 268, 268, 268, 268, 269, 269, 269, 269, 269, 270, 270, 270, 270, 270, 271, 271, 271, 271, 271, 272, 272, 272, 272, 272, 273, 273, 273, 273, 273, 274, 274, 274, 274, 274, 275, 275, 275, 275, 275, 276, 276, 276, 276, 276, 277, 277, 277, 277, 277, 278, 278, 278, 278, 278, 279, 279, 279, 279, 279, 280, 280, 280, 280, 280, 281, 281, 281, 281, 281, 282, 282, 282, 282, 282, 283, 283, 283, 283, 283, 284, 284, 284, 284, 284, 285, 285, 285, 285, 285, 286, 286, 286, 286, 286, 287, 287, 287, 287, 287, 288, 288, 288, 288, 288, 289, 289, 289, 289, 289, 290, 290, 290, 290, 290, 291, 291, 291, 291, 291, 292, 292, 292, 292, 292, 293, 293, 293, 293, 293, 294, 294, 294, 294, 294, 295, 295, 295, 295, 295, 296, 296, 296, 296, 296, 297, 297, 297, 297, 297, 298, 298, 298, 298, 298, 299, 299, 299, 299, 299, 300, 300, 300, 300, 300, 301, 301, 301, 301, 301, 302, 302, 302, 302, 302, 303, 303, 303, 303, 303, 304, 304, 304, 304, 304, 305, 305, 305, 305, 305, 306, 306, 306, 306, 306, 307, 307, 307, 307, 307, 308, 308, 308, 308, 308, 309, 309, 309, 309, 309, 310, 310, 310, 310, 310, 311, 311, 311, 311, 311, 312, 312, 312, 312, 312, 313, 313, 313, 313, 313, 314, 314, 314, 314, 314, 315, 315, 315, 315, 315, 316, 316, 316, 316, 316, 317, 317, 317, 317, 317, 318, 318, 318, 318, 318, 319, 319, 319, 319, 319, 320, 320, 320, 320, 320, 321, 321, 321, 321, 321, 322, 322, 322, 322, 322, 323, 323, 323, 323, 323, 324, 324, 324, 324, 324, 325, 325, 325, 325, 325, 326, 326, 326, 326, 326, 327, 327, 327, 327, 327, 328, 328, 328, 328, 328, 329, 329, 329, 329, 329, 330, 330, 330, 330, 330, 331, 331, 331, 331, 331, 332, 332, 332, 332, 332, 333, 333, 333, 333, 333, 334, 334, 334, 334, 334, 335, 335, 335, 335, 335, 336, 336, 336, 336, 336, 337, 337, 337, 337, 337, 338, 338, 338, 338, 338, 339, 339, 339, 339, 339, 340, 340, 340, 340, 340, 341, 341, 341, 341, 341, 342, 342, 342, 342, 342, 343, 343, 343, 343, 343, 344, 344, 344, 344, 344, 345, 345, 345, 345, 345, 346, 346, 346, 346, 346, 347, 347, 347, 347, 347, 348, 348, 348, 348, 348, 349, 349, 349, 349, 349, 350, 350, 350, 350, 350, 351, 351, 351, 351, 351, 352, 352, 352, 352, 352, 353, 353, 353, 353, 353, 354, 354, 354, 354, 354, 355, 355, 355, 355, 355, 356, 356, 356, 356, 356, 357, 357, 357, 357, 357, 358, 358, 358, 358, 358, 359, 359, 359, 359, 359, 360, 360, 360, 360, 360, 361, 361, 361, 361, 361, 362, 362, 362, 362, 362, 363, 363, 363, 363, 363, 364, 364, 364, 364, 364, 365, 365, 365, 365, 365, 366, 366, 366, 366, 366, 367, 367, 367, 367, 367, 368, 368, 368, 368, 368, 369, 369, 369, 369, 369, 370, 370, 370, 370, 370, 371, 371, 371, 371, 371, 372, 372, 372, 372, 372, 373, 373, 373, 373, 373, 374, 374, 374, 374, 374, 375, 375, 375, 375, 375, 376, 376, 376, 376, 376, 377, 377, 377, 377, 377, 378, 378, 378, 378, 378, 379, 379, 379, 379, 379, 380, 380, 380, 380, 380, 381, 381, 381, 381, 381, 382, 382, 382, 382, 382, 383, 383, 383, 383, 383, 384, 384, 384, 384, 384, 385, 385, 385, 385, 385, 386, 386, 386, 386, 386, 387, 387, 387, 387, 387, 388, 388, 388, 388, 388, 389, 389, 389, 389, 389, 390, 390, 390, 390, 390, 391, 391, 391, 391, 391, 392, 392, 392, 392, 392, 393, 393, 393, 393, 393, 394, 394, 394, 394, 394, 395, 395, 395, 395, 395, 396, 396, 396, 396, 396, 397, 397, 397, 397, 397, 398, 398, 398, 398, 398, 399, 399, 399, 399, 399, 400, 400, 400, 400, 400, 401, 401, 401, 401, 401, 402, 402, 402, 402, 402, 403, 403, 403, 403, 403, 404, 404, 404, 404, 404, 405, 405, 405, 405, 405, 406, 406, 406, 406, 406, 407, 407, 407, 407, 407, 408, 408, 408, 408, 408, 409, 409, 409, 409, 409, 410, 410, 410, 410, 410, 411, 411, 411, 411, 411, 412, 412, 412, 412, 412, 413, 413, 413, 413, 413, 414, 414, 414, 414, 414, 415, 415, 415, 415, 415, 416, 416, 416, 416, 416, 417, 417, 417, 417, 417, 418, 418, 418, 418, 418, 419, 419, 419, 419, 419, 420, 420, 420, 420, 420, 421, 421, 421, 421, 421, 422, 422, 422, 422, 422, 423, 423, 423, 423, 423, 424, 424, 424, 424, 424, 425, 425, 425, 425, 425, 426, 426, 426, 426, 426, 427, 427, 427, 427, 427, 428, 428, 428, 428, 428, 429, 429, 429, 429, 429, 430, 430, 430, 430, 430, 431, 431, 431, 431, 431, 432, 432, 432, 432, 432, 433, 433, 433, 433, 433, 434, 434, 434, 434, 434, 435, 435, 435, 435, 435, 436, 436, 436, 436, 436, 437, 437, 437, 437, 437, 438, 438, 438, 438, 438, 439, 439, 439, 439, 439, 440, 440, 440, 440, 440, 441, 441, 441, 441, 441, 442, 442, 442, 442, 442, 443, 443, 443, 443, 443, 444, 444, 444, 444, 444, 445, 445, 445, 445, 445, 446, 446, 446, 446, 446, 447, 447, 447, 447, 447, 448, 448, 448, 448, 448, 449, 449, 449, 449, 449, 450, 450, 450, 450, 450, 451, 451, 451, 451, 451, 452, 452, 452, 452, 452, 453, 453, 453, 453, 453, 454, 454, 454, 454, 454, 455, 455, 455, 455, 455, 456, 456, 456, 456, 456, 457, 457, 457, 457, 457, 458, 458, 458, 458, 458, 459, 459, 459, 459, 459, 460, 460, 460, 460, 460, 461, 461, 461, 461, 461, 462, 462, 462, 462, 462, 463, 463, 463, 463, 463, 464, 464, 464, 464, 464, 465, 465, 465, 465, 465, 466, 466, 466, 466, 466, 467, 467, 467, 467, 467, 468, 468, 468, 468, 468, 469, 469, 469, 469, 469, 470, 470, 470, 470, 470, 471, 471, 471, 471, 471, 472, 472, 472, 472, 472, 473, 473, 473, 473, 473, 474, 474, 474, 474, 474, 475, 475, 475, 475, 475, 476, 476, 476, 476, 476, 477, 477, 477, 477, 477, 478, 478, 478, 478, 478, 479, 479, 479, 479, 479, 480, 480, 480, 480, 480, 481, 481, 481, 481, 481, 482, 482, 482, 482, 482, 483, 483, 483, 483, 483, 484, 484, 484, 484, 484, 485, 485, 485, 485, 485, 486, 486, 486, 486, 486, 487, 487, 487, 487, 487, 488, 488, 488, 488, 488, 489, 489, 489, 489, 489, 490, 490, 490, 490, 490, 491, 491, 491, 491, 491, 492, 492, 492, 492, 492, 493, 493, 493, 493, 493, 494, 494, 494, 494, 494, 495, 495, 495, 495, 495, 496, 496, 496, 496, 496, 497, 497, 497, 497, 497, 498, 498, 498, 498, 498, 499, 499, 499, 499, 499, 500, 500, 500, 500, 500, 501, 501, 501, 501, 501, 502, 502, 502, 502, 502, 503, 503, 503, 503, 503, 504, 504, 504, 504, 504, 505, 505, 505, 505, 505, 506, 506, 506, 506, 506, 507, 507, 507, 507, 507, 508, 508, 508, 508, 508, 509, 509, 509, 509, 509, 510, 510, 510, 510, 510, 511, 511, 511, 511, 511, 512, 512, 512, 512, 512, 513, 513, 513, 513, 513, 514, 514, 514, 514, 514, 515, 515, 515, 515, 515, 516, 516, 516, 516, 516, 517, 517, 517, 517, 517, 518, 518, 518, 518, 518, 519, 519, 519, 519, 519, 520, 520, 520, 520, 520, 521, 521, 521, 521, 521, 522, 522, 522, 522, 522, 523, 523, 523, 523, 523, 524, 524, 524, 524, 524, 525, 525, 525, 525, 525, 526, 526, 526, 526, 526, 527, 527, 527, 527, 527, 528, 528, 528, 528, 528, 529, 529, 529, 529, 529, 530, 530, 530, 530, 530, 531, 531, 531, 531, 531, 532, 532, 532, 532, 532, 533, 533, 533, 533, 533, 534, 534, 534, 534, 534, 535, 535, 535, 535, 535, 536, 536, 536, 536, 536, 537, 537, 537, 537, 537, 538, 538, 538, 538, 538, 539, 539, 539, 539, 539, 540, 540, 540, 540, 540, 541, 541, 541, 541, 541, 542, 542, 542, 542, 542, 543, 543, 543, 543, 543, 544, 544, 544, 544, 544, 545, 545, 545, 545, 545, 546, 546, 546, 546, 546, 547, 547, 547, 547, 547, 548, 548, 548, 548, 548, 549, 549, 549, 549, 549, 550, 550, 550, 550, 550, 551, 551, 551, 551, 551, 552, 552, 552, 552, 552, 553, 553, 553, 553, 553, 554, 554, 554, 554, 554, 555, 555, 555, 555, 555, 556, 556, 556, 556, 556, 557, 557, 557, 557, 557, 558, 558, 558, 558, 558, 559, 559, 559, 559, 559, 560, 560, 560, 560, 560, 561, 561, 561, 561, 561, 562, 562, 562, 562, 562, 563, 563, 563, 563, 563, 564, 564, 564, 564, 564, 565, 565, 565, 565, 565, 566, 566, 566, 566, 566, 567, 567, 567, 567, 567, 568, 568, 568, 568, 568, 569, 569, 569, 569, 569, 570, 570, 570, 570, 570, 571, 571, 571, 571, 571, 572, 572, 572, 572, 572, 573, 573, 573, 573, 573, 574, 574, 574, 574, 574, 575, 575, 575, 575, 575, 576, 576, 576, 576, 576, 577, 577, 577, 577, 577, 578, 578, 578, 578, 578, 579, 579, 579, 579, 579, 580, 580, 580, 580, 580, 581, 581, 581, 581, 581, 582, 582, 582, 582, 582, 583, 583, 583, 583, 583, 584, 584, 584, 584, 584, 585, 585, 585, 585, 585, 586, 586, 586, 586, 586, 587, 587, 587, 587, 587, 588, 588, 588, 588, 588, 589, 589, 589, 589, 589, 590, 590, 590, 590, 590, 591, 591, 591, 591, 591, 592, 592, 592, 592, 592, 593, 593, 593, 593, 593, 594, 594, 594, 594, 594, 595, 595, 595, 595, 595, 596, 596, 596, 596, 596, 597, 597, 597, 597, 597, 598, 598, 598, 598, 598, 599, 599, 599, 599, 599, 600, 600, 600, 600, 600, 601, 601, 601, 601, 601, 602, 602, 602, 602, 602, 603, 603, 603, 603, 603, 604, 604, 604, 604, 604, 605, 605, 605, 605, 605, 606, 606, 606, 606, 606, 607, 607, 607, 607, 607, 608, 608, 608, 608, 608, 609, 609, 609, 609, 609, 610, 610, 610, 610, 610, 611, 611, 611, 611, 611, 612, 612, 612, 612, 612, 613, 613, 613, 613, 613, 614, 614, 614, 614, 614, 615, 615, 615, 615, 615, 616, 616, 616, 616, 616, 617, 617, 617, 617, 617, 618, 618, 618, 618, 618, 619, 619, 619, 619, 619, 620, 620, 620, 620, 620, 621, 621, 621, 621, 621, 622, 622, 622, 622, 622, 623, 623, 623, 623, 623, 624, 624, 624, 624, 624, 625, 625, 625, 625, 625, 626, 626, 626, 626, 626, 627, 627, 627, 627, 627, 628, 628, 628, 628, 628, 629, 629, 629, 629, 629, 630, 630, 630, 630, 630, 631, 631, 631, 631, 631, 632, 632, 632, 632, 632, 633, 633, 633, 633, 633, 634, 634, 634, 634, 634, 635, 635, 635, 635, 635, 636, 636, 636, 636, 636, 637, 637, 637, 637, 637, 638, 638, 638, 638, 638, 639, 639, 639, 639, 639, 640, 640, 640, 640, 640, 641, 641, 641, 641, 641, 642, 642, 642, 642, 642, 643, 643, 643, 643, 643, 644, 644, 644, 644, 644, 645, 645, 645, 645, 645, 646, 646, 646, 646, 646, 647, 647, 647, 647, 647, 648, 648, 648, 648, 648, 649, 649, 649, 649, 649, 650, 650, 650, 650, 650, 651, 651, 651, 651, 651, 652, 652, 652, 652, 652, 653, 653, 653, 653, 653, 654, 654, 654, 654, 654, 655, 655, 655, 655, 655, 656, 656, 656, 656, 656, 657, 657, 657, 657, 657, 658, 658, 658, 658, 658, 659, 659, 659, 659, 659, 660, 660, 660, 660, 660, 661, 661, 661, 661, 661, 662, 662, 662, 662, 662, 663, 663, 663, 663, 663, 664, 664, 664, 664, 664, 665, 665, 665, 665, 665, 666, 666, 666, 666, 666, 667, 667, 667, 667, 667, 668, 668, 668, 668, 668, 669, 669, 669, 669, 669, 670, 670, 670, 670, 670, 671, 671, 671, 671, 671, 672, 672, 672, 672, 672, 673, 673, 673, 673, 673, 674, 674, 674, 674, 674, 675, 675, 675, 675, 675, 676, 676, 676, 676, 676, 677, 677, 677, 677, 677, 678, 678, 678, 678, 678, 679, 679, 679, 679, 679, 680, 680, 680, 680, 680, 681, 681, 681, 681, 681, 682, 682, 682, 682, 682, 683, 683, 683, 683, 683, 684, 684, 684, 684, 684, 685, 685, 685, 685, 685, 686, 686, 686, 686, 686, 687, 687, 687, 687, 687, 688, 688, 688, 688, 688, 689, 689, 689, 689, 689, 690, 690, 690, 690, 690, 691, 691, 691, 691, 691, 692, 692, 692, 692, 692, 693, 693, 693, 693, 693, 694, 694, 694, 694, 694, 695, 695, 695, 695, 695, 696, 696, 696, 696, 696, 697, 697, 697, 697, 697, 698, 698, 698, 698, 698, 699, 699, 699, 699, 699, 700, 700, 700, 700, 700, 701, 701, 701, 701, 701, 702, 702, 702, 702, 702, 703, 703, 703, 703, 703, 704, 704, 704, 704, 704, 705, 705, 705, 705, 705, 706, 706, 706, 706, 706, 707, 707, 707, 707, 707, 708, 708, 708, 708, 708, 709, 709, 709, 709, 709, 710, 710, 710, 710, 710, 711, 711, 711, 711, 711, 712, 712, 712, 712, 712, 713, 713, 713, 713, 713, 714, 714, 714, 714, 714, 715, 715, 715, 715, 715, 716, 716, 716, 716, 716, 717, 717, 717, 717, 717, 718, 718, 718, 718, 718, 719, 719, 719, 719, 719, 720, 720, 720, 720, 720, 721, 721, 721, 721, 721, 722, 722, 722, 722, 722, 723, 723, 723, 723, 723, 724, 724, 724, 724, 724, 725, 725, 725, 725, 725, 726, 726, 726, 726, 726, 727, 727, 727, 727, 727, 728, 728, 728, 728, 728, 729, 729, 729, 729, 729, 730, 730, 730, 730, 730, 731, 731, 731, 731, 731, 732, 732, 732, 732, 732, 733, 733, 733, 733, 733, 734, 734, 734, 734, 734, 735, 735, 735, 735, 735, 736, 736, 736, 736, 736, 737, 737, 737, 737, 737, 738, 738, 738, 738, 738, 739, 739, 739, 739, 739, 740, 740, 740, 740, 740, 741, 741, 741, 741, 741, 742, 742, 742, 742, 742, 743, 743, 743, 743, 743, 744, 744, 744, 744, 744, 745, 745, 745, 745, 745, 746, 746, 746, 746, 746, 747, 747, 747, 747, 747, 748, 748, 748, 748, 748, 749, 749, 749, 749, 749, 750, 750, 750, 750, 750, 751, 751, 751, 751, 751, 752, 752, 752, 752, 752, 753, 753, 753, 753, 753, 754, 754, 754, 754, 754, 755, 755, 755, 755, 755, 756, 756, 756, 756, 756, 757, 757, 757, 757, 757, 758, 758, 758, 758, 758, 759, 759, 759, 759, 759, 760, 760, 760, 760, 760, 761, 761, 761, 761, 761, 762, 762, 762, 762, 762, 763, 763, 763, 763, 763, 764, 764, 764, 764, 764, 765, 765, 765, 765, 765, 766, 766, 766, 766, 766, 767, 767, 767, 767, 767, 768, 768, 768, 768, 768, 769, 769, 769, 769, 769, 770, 770, 770, 770, 770, 771, 771, 771, 771, 771, 772, 772, 772, 772, 772, 773, 773, 773, 773, 773, 774, 774, 774, 774, 774, 775, 775, 775, 775, 775, 776, 776, 776, 776, 776, 777, 777, 777, 777, 777, 778, 778, 778, 778, 778, 779, 779, 779, 779, 779, 780, 780, 780, 780, 780, 781, 781, 781, 781, 781, 782, 782, 782, 782, 782, 783, 783, 783, 783, 783, 784, 784, 784, 784, 784, 785, 785, 785, 785, 785, 786, 786, 786, 786, 786, 787, 787, 787, 787, 787, 788, 788, 788, 788, 788, 789, 789, 789, 789, 789, 790, 790, 790, 790, 790, 791, 791, 791, 791, 791, 792, 792, 792, 792, 792, 793, 793, 793, 793, 793, 794, 794, 794, 794, 794, 795, 795, 795, 795, 795, 796, 796, 796, 796, 796, 797, 797, 797, 797, 797, 798, 798, 798, 798, 798, 799, 799, 799, 799, 799, 800, 800, 800, 800, 800, 801, 801, 801, 801, 801, 802, 802, 802, 802, 802, 803, 803, 803, 803, 803, 804, 804, 804, 804, 804, 805, 805, 805, 805, 805, 806, 806, 806, 806, 806, 807, 807, 807, 807, 807, 808, 808, 808, 808, 808, 809, 809, 809, 809, 809, 810, 810, 810, 810, 810, 811, 811, 811, 811, 811, 812, 812, 812, 812, 812, 813, 813, 813, 813, 813, 814, 814, 814, 814, 814, 815, 815, 815, 815, 815, 816, 816, 816, 816, 816, 817, 817, 817, 817, 817, 818, 818, 818, 818, 818, 819, 819, 819, 819, 819, 820, 820, 820, 820, 820, 821, 821, 821, 821, 821, 822, 822, 822, 822, 822, 823, 823, 823, 823, 823, 824, 824, 824, 824, 824, 825, 825, 825, 825, 825, 826, 826, 826, 826, 826, 827, 827, 827, 827, 827, 828, 828, 828, 828, 828, 829, 829, 829, 829, 829, 830, 830, 830, 830, 830, 831, 831, 831, 831, 831, 832, 832, 832, 832, 832, 833, 833, 833, 833, 833, 834, 834, 834, 834, 834, 835, 835, 835, 835, 835, 836, 836, 836, 836, 836, 837, 837, 837, 837, 837, 838, 838, 838, 838, 838, 839, 839, 839, 839, 839, 840, 840, 840, 840, 840, 841, 841, 841, 841, 841, 842, 842, 842, 842, 842, 843, 843, 843, 843, 843, 844, 844, 844, 844, 844, 845, 845, 845, 845, 845, 846, 846, 846, 846, 846, 847, 847, 847, 847, 847, 848, 848, 848, 848, 848, 849, 849, 849, 849, 849, 850, 850, 850, 850, 850, 851, 851, 851, 851, 851, 852, 852, 852, 852, 852, 853, 853, 853, 853, 853, 854, 854, 854, 854, 854, 855, 855, 855, 855, 855, 856, 856, 856, 856, 856, 857, 857, 857, 857, 857, 858, 858, 858, 858, 858, 859, 859, 859, 859, 859, 860, 860, 860, 860, 860, 861, 861, 861, 861, 861, 862, 862, 862, 862, 862, 863, 863, 863, 863, 863, 864, 864, 864, 864, 864, 865, 865, 865, 865, 865, 866, 866, 866, 866, 866, 867, 867, 867, 867, 867, 868, 868, 868, 868, 868, 869, 869, 869, 869, 869, 870, 870, 870, 870, 870, 871, 871, 871, 871, 871, 872, 872, 872, 872, 872, 873, 873, 873, 873, 873, 874, 874, 874, 874, 874, 875, 875, 875, 875, 875, 876, 876, 876, 876, 876, 877, 877, 877, 877, 877, 878, 878, 878, 878, 878, 879, 879, 879, 879, 879, 880, 880, 880, 880, 880, 881, 881, 881, 881, 881, 882, 882, 882, 882, 882, 883, 883, 883, 883, 883, 884, 884, 884, 884, 884, 885, 885, 885, 885, 885, 886, 886, 886, 886, 886, 887, 887, 887, 887, 887, 888, 888, 888, 888, 888, 889, 889, 889, 889, 889, 890, 890, 890, 890, 890, 891, 891, 891, 891, 891, 892, 892, 892, 892, 892, 893, 893, 893, 893, 893, 894, 894, 894, 894, 894, 895, 895, 895, 895, 895, 896, 896, 896, 896, 896, 897, 897, 897, 897, 897, 898, 898, 898, 898, 898, 899, 899, 899, 899, 899, 900, 900, 900, 900, 900, 901, 901, 901, 901, 901, 902, 902, 902, 902, 902, 903, 903, 903, 903, 903, 904, 904, 904, 904, 904, 905, 905, 905, 905, 905, 906, 906, 906, 906, 906, 907, 907, 907, 907, 907, 908, 908, 908, 908, 908, 909, 909, 909, 909, 909, 910, 910, 910, 910, 910, 911, 911, 911, 911, 911, 912, 912, 912, 912, 912, 913, 913, 913, 913, 913, 914, 914, 914, 914, 914, 915, 915, 915, 915, 915, 916, 916, 916, 916, 916, 917, 917, 917, 917, 917, 918, 918, 918, 918, 918, 919, 919, 919, 919, 919, 920, 920, 920, 920, 920, 921, 921, 921, 921, 921, 922, 922, 922, 922, 922, 923, 923, 923, 923, 923, 924, 924, 924, 924, 924, 925, 925, 925, 925, 925, 926, 926, 926, 926, 926, 927, 927, 927, 927, 927, 928, 928, 928, 928, 928, 929, 929, 929, 929, 929, 930, 930, 930, 930, 930, 931, 931, 931, 931, 931, 932, 932, 932, 932, 932, 933, 933, 933, 933, 933, 934, 934, 934, 934, 934, 935, 935, 935, 935, 935, 936, 936, 936, 936, 936, 937, 937, 937, 937, 937, 938, 938, 938, 938, 938, 939, 939, 939, 939, 939, 940, 940, 940, 940, 940, 941, 941, 941, 941, 941, 942, 942, 942, 942, 942, 943, 943, 943, 943, 943, 944, 944, 944, 944, 944, 945, 945, 945, 945, 945, 946, 946, 946, 946, 946, 947, 947, 947, 947, 947, 948, 948, 948, 948, 948, 949, 949, 949, 949, 949, 950, 950, 950, 950, 950, 951, 951, 951, 951, 951, 952, 952, 952, 952, 952, 953, 953, 953, 953, 953, 954, 954, 954, 954, 954, 955, 955, 955, 955, 955, 956, 956, 956, 956, 956, 957, 957, 957, 957, 957, 958, 958, 958, 958, 958, 959, 959, 959, 959, 959, 960, 960, 960, 960, 960, 961, 961, 961, 961, 961, 962, 962, 962, 962, 962, 963, 963, 963, 963, 963, 964, 964, 964, 964, 964, 965, 965, 965, 965, 965, 966, 966, 966, 966, 966, 967, 967, 967, 967, 967, 968, 968, 968, 968, 968, 969, 969, 969, 969, 969, 970, 970, 970, 970, 970, 971, 971, 971, 971, 971, 972, 972, 972, 972, 972, 973, 973, 973, 973, 973, 974, 974, 974, 974, 974, 975, 975, 975, 975, 975, 976, 976, 976, 976, 976, 977, 977, 977, 977, 977, 978, 978, 978, 978, 978, 979, 979, 979, 979, 979, 980, 980, 980, 980, 980, 981, 981, 981, 981, 981, 982, 982, 982, 982, 982, 983, 983, 983, 983, 983, 984, 984, 984, 984, 984, 985, 985, 985, 985, 985, 986, 986, 986, 986, 986, 987, 987, 987, 987, 987, 988, 988, 988, 988, 988, 989, 989, 989, 989, 989, 990, 990, 990, 990, 990, 991, 991, 991, 991, 991, 992, 992, 992, 992, 992, 993, 993, 993, 993, 993, 994, 994, 994, 994, 994, 995, 995, 995, 995, 995, 996, 996, 996, 996, 996, 997, 997, 997, 997, 997, 998, 998, 998, 998, 998, 999, 999, 999, 999, 999, 1000, 1000, 1000, 1000, 1000, 1001, 1001, 1001, 1001, 1001, 1002, 1002, 1002, 1002, 1002, 1003, 1003, 1003, 1003, 1003, 1004, 1004, 1004, 1004, 1004, 1005, 1005, 1005, 1005, 1005, 1006, 1006, 1006, 1006, 1006, 1007, 1007, 1007, 1007, 1007, 1008, 1008, 1008, 1008, 1008, 1009, 1009, 1009, 1009, 1009, 1010, 1010, 1010, 1010, 1010, 1011, 1011, 1011, 1011, 1011, 1012, 1012, 1012, 1012, 1012, 1013, 1013, 1013, 1013, 1013, 1014, 1014, 1014, 1014, 1014, 1015, 1015, 1015, 1015, 1015, 1016, 1016, 1016, 1016, 1016, 1017, 1017, 1017, 1017, 1017, 1018, 1018, 1018, 1018, 1018, 1019, 1019, 1019, 1019, 1019, 1020, 1020, 1020, 1020, 1020, 1021, 1021, 1021, 1021, 1021, 1022, 1022, 1022, 1022, 1022, 1023, 1023, 1023, 1023, 1023, 1024, 1024, 1024, 1024, 1024], [4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0  …  -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0], 1024, 1024), [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0], [0.0 0.03842943919353914 … 0.15224093497742697 0.03842943919353936; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785; … ; 0.15224093497742697 0.1906703741709661 … 0.30448186995485393 0.19067037417096633; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785])

One of the benefits of the Bayesian approach is that we can fit for the hyperparameters of our prior/regularizers unlike traditional RML appraoches. To construct this heirarchical prior we will first make a map that takes in our regularizer hyperparameters and returns the image prior given those hyperparameters.

fmap = let crcache=crcache
    x->GaussMarkovRandomField(x, 1.0, crcache)
end
#1 (generic function with 1 method)

Now we can finally form our image prior. For this we use a heirarchical prior where the inverse correlation length is given by a Half-Normal distribution whose peak is at zero and standard deviation is 0.1/rat where recall rat is the beam size per pixel. For the variance of the random field we use another half normal prior with standard deviation 0.1. The reason we use the half-normal priors is to prefer "simple" structures. Gaussian Markov random fields are extremly flexible models, and to prevent overfitting it is common to use priors that penalize complexity. Therefore, we want to use priors that enforce similarity to our mean image. If the data wants more complexity then it will drive us away from the prior.

cprior = HierarchicalPrior(fmap, InverseGamma(1.0, -log(0.01*rat)))
VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)

)

We can now form our model parameter priors. Like our other imaging examples, we use a Dirichlet prior for our image pixels. For the log gain amplitudes, we use the CalPrior which automatically constructs the prior for the given jones cache gcache.

prior = NamedDist(
         c = cprior,
         fg = Uniform(0.0, 1.0),
         σimg = truncated(Normal(0.0, 0.1); lower = 0.0),
         lgamp = CalPrior(distamp, gcache),
         gphase = CalPrior(distphase, gcachep),
        )
VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)

)
, Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0), CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)

jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778  …  0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)

jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)

Putting it all together we form our likelihood and posterior objects for optimization and sampling.

lklhd = RadioLikelihood(sky, instrument, dvis; skymeta=metadata, instrumentmeta=metadata)
post = Posterior(lklhd, prior)
Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Comrade.ModelMetadata{typeof(Main.instrument), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#26#27"{Float64, Vector{Float64}}, Vector{ComplexF64}}}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}}(RadioLikelihood
	Number of data products: 1
, VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)

)
, Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0), CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)

jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778  …  0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)

jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)
)

Reconstructing the Image and Instrument Effects

To sample from this posterior, it is convenient to move from our constrained parameter space to an unconstrained one (i.e., the support of the transformed posterior is (-∞, ∞)). This is done using the asflat function.

tpost = asflat(post)
ndim = dimension(tpost)
1364

Our Posterior and TransformedPosterior objects satisfy the LogDensityProblems interface. This allows us to easily switch between different AD backends and many of Julia's statistical inference packages use this interface as well.

using LogDensityProblemsAD
using Zygote
gtpost = ADgradient(Val(:Zygote), tpost)
x0 = randn(rng, ndim)
LogDensityProblemsAD.logdensity_and_gradient(gtpost, x0)
(-1.0889677227706835e8, [3.971486365346438, -41.82940853200138, 69.59557261663501, 4.56980536621942, 10.789295400079078, 240.91028613995184, 39.14697994202466, 42.21429448053749, 50.17133106361702, 69.99422943458495  …  3666.7681059308334, 11816.75911472515, 1587.1490272240167, -343.3074914007028, 6760.359870612455, 349.8400750594293, 17534.646730431738, -7352.6390765479655, -181.37189823832506, 410.6393309543057])

We can now also find the dimension of our posterior or the number of parameters we are going to sample.

Warning

This can often be different from what you would expect. This is especially true when using angular variables where we often artificially increase the dimension of the parameter space to make sampling easier.

To initialize our sampler we will use optimize using LBFGS

using ComradeOptimization
using OptimizationOptimJL
f = OptimizationFunction(tpost, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(f, prior_sample(rng, tpost), nothing)
ℓ = logdensityof(tpost)
sol = solve(prob, LBFGS(), maxiters=4_000, g_tol=1e-1);

Now transform back to parameter space

xopt = transform(tpost, sol.u)
(c = (params = [-2.1098589629146646e-6 -2.7793645267644962e-6 … -1.2269581365980504e-5 -6.085786597700291e-7; -1.572557001092337e-5 -1.15714428165114e-5 … -1.3313614918637788e-5 3.6943445068705693e-6; … ; 5.079494316638737e-6 6.6231224439995475e-6 … 6.936664448300607e-6 3.705488339991044e-6; 7.696337466030705e-6 4.936412393809447e-6 … 8.662204871094582e-6 -2.969981442452276e-6], hyperparams = 0.05621916462568886), fg = 0.004543326958453666, σimg = 3.1097487638021843, lgamp = [-0.008547781982182092, 0.008536379590309503, -0.8869750654212027, 0.13574789624857495, 0.014653864801009768, -0.014176764899803103, -0.315855393599793, 0.24197386871827803, 0.00545369852654422, -0.005059738697923953  …  -0.29459460575971685, -0.0006657318984340658, -1.0921314321342268, -0.0029435606343759203, 0.02791213878701203, -0.02725007707447156, -0.2844298383950318, 0.000527583564805484, -1.1510095985232354, -0.004276382420868597], gphase = [-0.8142830755077827, -3.1087022029209685, 2.3940391487654473, -0.8711512360966408, -3.011060247261388, 2.200582319770767, -0.9285512324750242, -2.9246188743515775, 2.011587727945264, -0.9773339783202181  …  -0.5846187216263065, -0.4677141980974109, -0.49304252468741316, -0.668860788126376, 1.2783702979811462, -0.5238627407076003, -0.5805848385889851, -0.46325486145882727, -0.8818178220857682, 1.3257595666627136])
Warning

Fitting gains tends to be very difficult, meaning that optimization can take a lot longer. The upside is that we usually get nicer images.

First we will evaluate our fit by plotting the residuals

using Plots
residual(vlbimodel(post, xopt), dvis)
Example block output

These look reasonable, although there may be some minor overfitting. This could be improved in a few ways, but that is beyond the goal of this quick tutorial. Plotting the image, we see that we have a much cleaner version of the closure-only image from Imaging a Black Hole using only Closure Quantities.

import WGLMakie as CM
img = intensitymap(skymodel(post, xopt), fovx, fovy, 128, 128)
CM.image(img, axis=(xreversed=true, aspect=1, title="MAP Image"), colormap=:afmhot)

Because we also fit the instrument model, we can inspect their parameters. To do this, Comrade provides a caltable function that converts the flattened gain parameters to a tabular format based on the time and its segmentation.

gt = Comrade.caltable(gcachep, xopt.gphase)
plot(gt, layout=(3,3), size=(600,500))
Example block output

The gain phases are pretty random, although much of this is due to us picking a random reference station for each scan.

Moving onto the gain amplitudes, we see that most of the gain variation is within 10% as expected except LMT, which has massive variations.

gt = Comrade.caltable(gcache, exp.(xopt.lgamp))
plot(gt, layout=(3,3), size=(600,500))
Example block output

To sample from the posterior, we will use HMC, specifically the NUTS algorithm. For information about NUTS, see Michael Betancourt's notes.

Note

For our metric, we use a diagonal matrix due to easier tuning

However, due to the need to sample a large number of gain parameters, constructing the posterior is rather time-consuming. Therefore, for this tutorial, we will only do a quick preliminary run, and any posterior inferences should be appropriately skeptical.

using ComradeAHMC
metric = DiagEuclideanMetric(ndim)
chain, stats = sample(rng, post, AHMC(;metric, autodiff=Val(:Zygote)), 700; nadapts=500, init_params=xopt)
(NamedTuple{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{NamedTuple{(:params, :hyperparams), Tuple{Matrix{Float64}, Float64}}, Float64, Float64, Vector{Float64}, Vector{Float64}}}[(c = (params = [-0.0015725897927871553 0.007451302908807232 … 0.01169677440194074 -0.007004156508935446; 0.007158389998256451 0.013374052201153775 … 0.00015583830133340825 -0.0032954785299928006; … ; -0.005647836760377006 0.004441779600554347 … -0.004533034056422414 0.0016256328196930152; -0.009678105412239266 0.010534205713748279 … -0.0035458525174033384 -0.0033108682246471978], hyperparams = 0.05629064524965654), fg = 0.004594316773724774, σimg = 3.122843508904704, lgamp = [-0.009277744786952095, 0.008698267474007053, -0.8908133939255265, 0.13354652588998334, 0.021262009985193994, -0.021153654231317136, -0.32413623044782475, 0.22930459702267203, 0.004270804936020348, -0.003785567148153178  …  -0.29592022744403257, 0.0018989328821043972, -1.0993121811316917, 0.00167018942556604, 0.025577942970800627, -0.025091973046137438, -0.2863943128035272, -0.006929448654620688, -1.1549519723252852, 0.0032493453006481343], gphase = [-0.8173179084060926, -3.0923549615652446, 2.402217579209197, -0.8685743481954618, -3.0227387860628196, 2.209568194293659, -0.9280484885449602, -2.9261935014323353, 2.0150550660914233, -0.9787525402303876  …  -0.5869710271332244, -0.47421813114509626, -0.47965748859380225, -0.6690678679943606, 1.2892863075849286, -0.5219048266587085, -0.578262219599456, -0.4647003632497338, -0.8844685701655433, 1.3201554680529122]), (c = (params = [-0.0015725897927871553 0.007451302908807232 … 0.01169677440194074 -0.007004156508935446; 0.007158389998256451 0.013374052201153775 … 0.00015583830133340825 -0.0032954785299928006; … ; -0.005647836760377006 0.004441779600554347 … -0.004533034056422414 0.0016256328196930152; -0.009678105412239266 0.010534205713748279 … -0.0035458525174033384 -0.0033108682246471978], hyperparams = 0.05629064524965654), fg = 0.004594316773724774, σimg = 3.122843508904704, lgamp = [-0.009277744786952095, 0.008698267474007053, -0.8908133939255265, 0.13354652588998334, 0.021262009985193994, -0.021153654231317136, -0.32413623044782475, 0.22930459702267203, 0.004270804936020348, -0.003785567148153178  …  -0.29592022744403257, 0.0018989328821043972, -1.0993121811316917, 0.00167018942556604, 0.025577942970800627, -0.025091973046137438, -0.2863943128035272, -0.006929448654620688, -1.1549519723252852, 0.0032493453006481343], gphase = [-0.8173179084060926, -3.0923549615652446, 2.402217579209197, -0.8685743481954618, -3.0227387860628196, 2.209568194293659, -0.9280484885449602, -2.9261935014323353, 2.0150550660914233, -0.9787525402303876  …  -0.5869710271332244, -0.47421813114509626, -0.47965748859380225, -0.6690678679943606, 1.2892863075849286, -0.5219048266587085, -0.578262219599456, -0.4647003632497338, -0.8844685701655433, 1.3201554680529122]), (c = (params = [-0.11452131855247688 -0.01382436452443098 … -0.1701782956337741 -0.1292608050929214; 0.012512144227187612 -0.024007508596418137 … -0.04997575931204093 -0.02397045499511187; … ; -0.08887760627023884 0.0010538975191178377 … -0.1051251312414246 -0.09346987252922477; -0.10911811081292389 -0.03564820235045443 … -0.05830644535718606 0.06186687718291096], hyperparams = 0.11935591939043778), fg = 0.003677418405085104, σimg = 2.545634516737658, lgamp = [0.0018671759340980112, -0.004673879892755381, -0.8821521715145803, 0.11233354276636909, -0.0013248077094337011, 0.003151674729635992, -0.31023214700223883, 0.24936839130230534, 0.0004408410412846085, -0.0013385537136734373  …  -0.3019551795325553, -0.0064092987928032995, -1.1114911159959497, -0.01119095669145198, 0.027528233628613515, -0.026306082803104154, -0.27910940032833964, -0.0006441452459057448, -1.1619795561083177, -0.003099218636952883], gphase = [-0.81871022335209, -3.127068335989979, 2.4194313098151845, -0.8665225912285376, -3.008587294561622, 2.224785801581465, -0.9292956506525154, -2.9218477299115504, 2.0360221209079845, -0.9776032963287072  …  -0.5811365240337019, -0.4753553589123778, -0.5153921769308379, -0.6709283150392051, 1.2561137644951412, -0.5211157153716282, -0.5871943952829662, -0.4976638838476991, -0.8764271545727449, 1.282727230599266]), (c = (params = [0.04276243924408368 0.09737196226667519 … -0.3057414673069623 0.0019055116521954148; -0.025518645685684026 -0.19090971659564873 … -0.16684701904070973 0.050086950346821514; … ; 0.032286920340220765 -0.04070765953927833 … -0.06549359807981041 -0.1096169937528537; -0.14566293124807067 -0.09085218503578937 … 0.010862799643638897 0.11680408197649061], hyperparams = 0.16934762761884778), fg = 0.004517975627366718, σimg = 2.019001706828473, lgamp = [0.0018921754245193837, -0.0015525108688682861, -0.9896351377289647, 0.0700422708917987, -0.006593109130283133, 0.002094300489605947, -0.3842520872175946, 0.20499439147749196, -0.008870138805901633, 0.003159672346234357  …  -0.2926053546694022, -0.017816728432547786, -1.1126896187123902, 0.0013330950815489747, 0.009898574203964962, -0.015001653151069454, -0.2706672166352087, 0.008285797467933984, -1.1665599829894135, -0.010615074443789528], gphase = [-0.8145583941268367, -3.062365385781824, 2.4101200514473535, -0.867465060275225, -2.979012187232555, 2.2212546570786182, -0.9335841237590328, -2.91879492281987, 2.03393203763219, -0.977333115888234  …  -0.587407869516936, -0.45767076697717873, -0.5571109816241582, -0.665343212969669, 1.2146892567536922, -0.5273802960137455, -0.5833943279391907, -0.500149354402202, -0.8773544937929335, 1.2832783622968391]), (c = (params = [-0.1206634440313447 -0.13379943951479598 … 0.35949114301434026 0.04987574834100288; 0.04045230964048701 0.22583136149721147 … 0.15495137277715773 -0.11631481532513957; … ; 0.026716566975185287 -0.007564283711263868 … 0.015253375488106671 0.16252207856433293; 0.23060544130597568 0.07547122006284036 … -0.07912227086834182 -0.08223473638529302], hyperparams = 0.24719367854112484), fg = 0.015106604581883946, σimg = 1.5771137000634798, lgamp = [0.009511593154449712, -0.009765927756296174, -1.1022582590381764, -0.0467715872024029, -0.0028249990124665637, 0.004257024675893884, -0.4995598966413009, 0.10895557220941758, -0.01867638456787596, 0.01931561880810533  …  -0.2689984479357474, 0.011557816770423508, -1.1359497011314106, -0.01422725893260789, -0.004386792090281303, 0.003082228857455045, -0.25162549585357513, -0.005891847199013752, -1.2043983597340011, 0.004975416035629001], gphase = [-0.813812775628389, -3.104722274162015, 2.2725200949396465, -0.871126308292479, -2.9954746811437567, 2.0693539140286434, -0.9309658011749881, -2.9331456789227235, 1.8732532540773368, -0.973643601734974  …  -0.5809614166862314, -0.478327361626811, -0.39880880490142556, -0.664491849448669, 1.367983322175391, -0.5305596109674162, -0.5905228249953516, -0.3407100065674994, -0.8757946511843301, 1.4462526602059442]), (c = (params = [0.04941551639361871 0.1315307974681064 … -0.4405343466404579 -0.12156325411305371; -0.05992036701396353 -0.3613122012192081 … -0.19638275061692714 0.05077699593506902; … ; -0.05402267077101361 -0.00894748080865622 … -0.037755796488889784 -0.20530725982750195; -0.2750668533977829 -0.11941996701063903 … 0.07277716618765871 0.056886664270239966], hyperparams = 0.2532317807975833), fg = 0.03218342355061539, σimg = 1.4293513378257934, lgamp = [-0.03290413989864965, 0.033580235773725015, -1.0677335668890748, 0.00822188849745387, -0.0019155590192328738, 0.00376303655389257, -0.48082027098960556, 0.12285903292210011, -0.019044543858051367, 0.021541743828785178  …  -0.2527089204680411, -0.014990975868634798, -1.130741079365554, 0.005974111565987598, -0.011903588594262852, 0.010174998353061859, -0.22313518749924832, 0.014616091468374405, -1.1710043972395732, -0.0015519903560844736], gphase = [-0.8068750965924057, -3.010613800370075, 2.3648226604456806, -0.8759876963116987, -2.9731998383562606, 2.1868616402609056, -0.9269513442178607, -2.915979673076525, 2.0030504873263837, -0.9787837656006282  …  -0.5872711194726661, -0.48788324159182134, -0.47435199399163513, -0.6606978160042588, 1.3127838060355779, -0.5213722142072325, -0.5853493888409413, -0.44081406024157943, -0.8913369848081604, 1.3378139992856948]), (c = (params = [0.03706165609885795 -0.30954175918519256 … -0.12930465294849058 0.38029400580951866; -0.2677666191509689 -0.5778374796090868 … -0.292514383600808 0.04151027613707038; … ; 0.13755686734661748 -0.4812933897109464 … -0.15733917501143813 -0.11000409487662816; 0.02104539930917976 -0.3237644787014109 … 0.016132435117229033 0.5970143295990613], hyperparams = 0.3308345550610656), fg = 0.034956200655189416, σimg = 1.6361777532950237, lgamp = [-0.02539311322198802, 0.024032142897012113, -1.0907079786165212, 0.06294949748316131, -0.0350107088323914, 0.02976291043164944, -0.47523979567979585, 0.2011634658705269, -0.021983083603930084, 0.02778464435904803  …  -0.25766396574191874, -0.013513197783121846, -1.1198645853178468, -0.0026400007000057057, -0.010278744349927168, 0.009040033706095345, -0.2432982451609293, 0.016674820387242016, -1.1804236001294626, -0.0006794470626027983], gphase = [-0.8065755769392506, -2.9956615890852802, 2.30060048006814, -0.8677601870925977, -2.9732135553007026, 2.1221627828027216, -0.9283051185894295, -2.9106688744058413, 1.9348537343808139, -0.975007003219888  …  -0.5814684348928316, -0.48155844596268677, -0.4069400084451019, -0.6546085092582716, 1.3322296106942677, -0.5275915830727714, -0.584680142724091, -0.3986108330439186, -0.8475516114372892, 1.4023251750215455]), (c = (params = [0.43022892822704534 -0.24852421474130001 … 0.6018835135379746 -0.16009489972961236; 0.133440021872474 0.11886763930615289 … 0.23157657495076514 -0.071542126469204; … ; -0.7840986034278833 0.0698699883612847 … -0.1249929837470144 -0.12509864921171218; -0.2876187619918358 -0.10740166699690705 … -0.37164363362898467 -0.5436167398696131], hyperparams = 0.4683746425187825), fg = 0.053547181944303975, σimg = 1.294803928675915, lgamp = [-0.01602247366875804, 0.012689862079915664, -1.1196196381741, -0.008968779559928823, -0.026196312077300693, 0.02637050340494914, -0.48340161085678957, 0.12363558170919788, -0.033722648438884235, 0.03077635656922559  …  -0.2447966872646166, 0.0026661704648207113, -1.102576873603745, -0.00768481456110504, -0.005338660385541464, 0.004128880374931566, -0.23702222075454552, 0.015091864585151353, -1.165597583998405, -0.019997605774312758], gphase = [-0.8071275988586573, -3.000476483792019, 2.1311775169104368, -0.8707709492915244, -2.916784660644283, 1.9619331177252066, -0.9293932725701534, -2.852436207878833, 1.762526241794694, -0.9806231681508146  …  -0.5886130845035851, -0.43264262339272236, -0.28019531624955324, -0.6400626452545839, 1.4799453613461027, -0.5249729290753788, -0.5374800365265999, -0.25937052122118015, -0.8638061397021286, 1.524760154901893]), (c = (params = [-0.21602927767905053 0.353113236537776 … 0.947992618973093 0.7490323956433578; -0.17934632450754406 -0.017369611956066133 … -0.4553386312894916 0.31684400221741565; … ; -0.5285176448941512 0.14389314792675062 … 0.6014413662122408 0.08885055759948339; 0.047250294493335405 0.23037192793568498 … 0.10977217968963728 0.7434801908015022], hyperparams = 0.5420476463242968), fg = 0.06516402927914762, σimg = 1.3619327442809042, lgamp = [-0.013893664623563588, 0.0056203856318408725, -1.1453559582032873, -0.005692369686719604, -0.01833876195919511, 0.016364511833559955, -0.52787118706063, 0.12320920199678587, -0.03445886980066742, 0.03716965951079681  …  -0.2570690282416446, 2.31918007965791e-5, -1.0894043193149077, -0.005194577028038992, -0.0007544603713095059, -0.0030758986333111657, -0.24757736163862304, 0.015056065597429337, -1.1561416726703524, -0.017866833135017333], gphase = [-0.8169822968349993, -3.012654325528695, 2.2236990968409076, -0.8718743667483115, -2.94334558437182, 2.0576506097232468, -0.9282877511432898, -2.8776927568249264, 1.8590310093684301, -0.9773695144888023  …  -0.5814793412169915, -0.44779424443568694, -0.35797929883308094, -0.6256830104429565, 1.4169431043830052, -0.5238983527408582, -0.557175662723935, -0.3243370975685813, -0.8351111433141247, 1.4700356434145534]), (c = (params = [0.3010374819923765 -0.47327557424829053 … -0.8316056837961368 -0.46485449326920314; -0.09934363531837198 -0.5056646155447402 … 0.9678112018429672 -0.04419019300101377; … ; 0.5585510635147889 -0.07446428964736293 … -0.8356261121413061 -0.3397323371592053; 0.046865478383529356 0.0006039702063531571 … -0.09245775836194678 -0.8312919226718025], hyperparams = 1.2095429856254747), fg = 0.46668150141627546, σimg = 1.1568236398815603, lgamp = [-0.043823763514277136, 0.03966270336528503, -0.8205822532438505, 0.03084826752205007, -0.003625556046066576, 0.0005057383860194079, -0.249868701980548, 0.128731718639066, -0.0446446753938895, 0.045458002619869584  …  -0.0280931610230094, -0.00845308526708314, -0.8322365791365242, -0.005017601281317137, 0.0021877792533397293, -0.0007796560848664699, -0.013842653356287742, 0.016222885508998895, -0.890432137180502, -0.025608980019298067], gphase = [-0.8159586819787127, -2.878704821153753, 2.3753269952819642, -0.8722803937919815, -2.835807306184451, 2.1806580846024186, -0.9291667000838109, -2.767265292497924, 1.989491384984162, -0.978078981494319  …  -0.5867990392916556, -0.32587015032352556, -0.31366095330152505, -0.5707687452308584, 1.4605379391998718, -0.5263364103319212, -0.4467676666044542, -0.2555477003338578, -0.7562956934786935, 1.5240407645323955])  …  (c = (params = [-0.8671370971444401 -0.3407681503715883 … -0.43051298642046365 -1.028481910683817; -1.8062668459663085 -1.0103230251357038 … -0.542203529780662 -1.4090156413971002; … ; -0.04585000344013028 -1.551780603990974 … -0.22073127723117852 0.14211322750448538; 0.11702669050174402 -0.9599309904424982 … 0.1680749879898696 -0.404749965893147], hyperparams = 12.879631145353983), fg = 0.35060672723073966, σimg = 1.065481960893482, lgamp = [-0.030257583724141317, 0.021104107701466666, -0.9154258298381848, 0.015246013345560652, 0.0025298432296797383, -0.0016350275576982266, -0.3600304739776619, 0.10825045131109506, -0.0323114795703467, 0.03257073225841522  …  -0.057939042812304796, 0.0020921017396304965, -0.907083500536196, -0.002611365459263511, -0.023700650689959656, 0.024425994281060136, -0.07879055884305088, 0.0058440821804296335, -0.9865061246925859, -0.005770029930313497], gphase = [-0.8059768804196814, 2.810181219886886, 0.3896932356713949, -0.8693080500440922, 2.9172606252270548, 0.16598685662563004, -0.9343517420075567, 2.957803537931578, -0.050004519562642395, -0.9783133520154312  …  -0.5872243787615846, -0.8948392015062353, 0.40872951420847026, -1.0210067015553732, 2.1849176922638662, -0.517791971083518, -0.9986183781471261, 0.48519996941224736, -1.2064640016477828, 2.2697236939373076]), (c = (params = [-0.02052417392752575 -0.415899014850036 … 0.06855073788047872 0.8583832461028875; -0.10509159312624411 -0.188203642728025 … -0.31269219228366485 0.06885957718607652; … ; -0.4629190613012997 0.46728089006936485 … -0.263798732245485 -0.003939485029152786; 0.11430075591894422 0.25604858099203764 … 0.4506988879850188 -0.12859574509855975], hyperparams = 14.393921130271638), fg = 0.3203954990755933, σimg = 1.022757451096168, lgamp = [-0.024269488680064958, 0.029335183425384388, -0.9218076879471147, -0.03158857682238199, -0.009121208715669177, 0.015215636900852678, -0.32705265446720766, 0.09739602511583559, -0.01732388007414719, 0.01583754360228864  …  -0.10784427966651335, -0.00030842332882258813, -0.9328831375664951, -0.023800894285911037, -0.01641859423041416, 0.014379415327771086, -0.10263716183361615, -0.0018243922363338505, -0.9651323459424452, 0.0009779307287550614], gphase = [-0.8162588371690299, 2.924198702701188, 0.6425245840511029, -0.8716024838157465, 2.988175076221865, 0.42549597886904067, -0.9247779380633041, 3.0360898591502785, 0.21158886324479467, -0.9763606572554538  …  -0.585534472123696, -0.8118065009300331, 0.2701545345539151, -0.988027982533267, 2.0509405190673933, -0.529076138406656, -0.9561699841539822, 0.31661505949070834, -1.1944304251237252, 2.1009680603001835]), (c = (params = [0.802791025465772 0.5714698729723124 … -0.2157700475497656 0.350369635690962; -0.3271130480895776 -0.64102778434859 … -0.621744145425193 0.23591280956742428; … ; 0.17614002974641071 -0.2862858101500911 … -0.04001031351192706 -0.06500627619065066; -0.14021328743139766 -0.09989428390988506 … -1.3644303645975862 0.19700616845954622], hyperparams = 27.95769517961656), fg = 0.3264645680316062, σimg = 0.9801806769320057, lgamp = [-0.03265684877697828, 0.03229080378557571, -0.8864303988279773, -0.002979056871012949, 0.004242435489936653, -0.01046551178945138, -0.3748246573665435, 0.09248111546888792, 0.0063603196903058426, -0.004217650721723957  …  -0.09409589063507963, 0.009029764641778489, -0.9023366708902171, -0.025531846774651795, -0.023004891086942106, 0.01960052858012735, -0.08308111593853494, -0.0036595481339312067, -0.9340124224296964, -0.0008542664726418009], gphase = [-0.8142951289432622, 2.9198026858701125, 0.6845160684975282, -0.8753001214947371, 2.985099221328402, 0.47113140377355006, -0.9326356829691207, 3.0317097093040832, 0.2698036092856532, -0.979020644793145  …  -0.5793177315322744, -0.7864810821314038, 0.26056742284286943, -0.9418542196518087, 2.0264101820676768, -0.5220342113603853, -0.9329673032474791, 0.31019502942033506, -1.139254684334117, 2.086061208679693]), (c = (params = [-0.17764053831079937 -0.4613152859359276 … 0.030098560378516374 0.19648447209259146; 0.056635408542460255 -0.5843678635676482 … -0.32402715793873194 -0.2506104076248472; … ; -0.6051724761489667 0.2075586479816824 … -0.3216992966987686 -0.4215484286742988; -0.12353050375187541 -0.15753504332148605 … 0.6818175381414371 -0.5423144418966654], hyperparams = 14.30507372862622), fg = 0.3407896112349495, σimg = 0.9541681142935324, lgamp = [-0.037660982157069964, 0.037786322289296, -0.9335528506498085, -0.00255118794136189, -0.0025027833955054406, 0.004428969157803507, -0.35459811584881706, 0.11975156035350105, -0.028466447177720636, 0.02382298303283026  …  -0.10411847625215671, 0.007269926782706109, -0.936401074332783, -0.012756603346552929, -0.01552768515003585, 0.02465265227416139, -0.08740187831786946, -0.0030421882649432283, -0.9700891117506597, 0.004549434834385142], gphase = [-0.8097002189947196, 2.895915587509408, 0.7910145601089971, -0.8751273920110642, 2.9721845529771898, 0.5694095371269537, -0.9352404288698906, 3.0179365831901657, 0.3444122446215692, -0.9738959379076301  …  -0.5860673620422845, -0.8007388832173545, 0.15557979360708474, -0.9040632459427682, 1.9278620397447876, -0.5224807625774129, -0.9253132651016879, 0.18884612431265468, -1.120939673662947, 1.9813934823260646]), (c = (params = [0.06322677364100741 -0.756068683478213 … -0.01864937524631345 0.10482321616896925; -0.48833527396632836 -0.33474724718360427 … -0.7745749643410267 -0.6314866927082876; … ; 0.4920341905848537 -0.4784191569344422 … -0.6851307640866121 -0.47231241376388094; 0.25679831889376203 -0.6611385834014166 … -1.0232592959448856 0.573169478101962], hyperparams = 13.263476004023902), fg = 0.3232546510172162, σimg = 0.9285679708673833, lgamp = [-0.008647481295217935, 0.008912963784199673, -0.8544456847030317, -0.002525482532072001, 0.005069479809285109, -0.0006063025448772048, -0.32614990038075165, 0.10799107092094229, -0.012200541471936093, 0.012939555730265229  …  -0.06487039518326854, -0.01469584861775613, -0.9457817418809483, -0.0035498159208012165, -0.020944520402137963, 0.022080120615868983, -0.0911969994710066, -0.0008797660375910779, -0.9896408438185932, 0.001788402172481854], gphase = [-0.8194734804758638, 2.8167349037075127, 0.5936592717918295, -0.8679100974868224, 2.8882163437139687, 0.3917048528979011, -0.9244684524851626, 2.9605181225745527, 0.16515314218490224, -0.978972780742266  …  -0.5837707442054967, -0.8189457195371035, 0.17090217598682678, -0.943405541988855, 1.9492413910794293, -0.5209000937338691, -0.9526956888872825, 0.22542600678164856, -1.1521297560627344, 2.0127164432028906]), (c = (params = [-0.15162579792051972 -0.08535972241732095 … -1.1978093079290368 -0.6778060435666774; 0.7583009791211905 0.02828378467991366 … -0.8066597321518204 0.03417282607875989; … ; -0.9805007587769419 -0.2530375585784624 … -0.21935719615099938 -0.16029780982485511; -0.8379566749704973 -0.3574051222516877 … 0.021655139459619725 -1.0675510164473847], hyperparams = 13.609326971519906), fg = 0.3244972413386148, σimg = 0.9948938559777447, lgamp = [-0.014315566295352144, 0.012785593155065457, -0.9167983944598221, 0.008147047614017343, -0.005496074355098304, 0.004532038374173607, -0.3639664709571025, 0.12688467961768965, -0.04311017393380548, 0.04988659366378765  …  -0.09512135076766308, 0.011726072321273873, -0.9629296970850436, -0.013810372653897361, 0.006617065122585911, -0.0034412944667943546, -0.11499063116893335, -0.0045277210837044906, -1.0187847707360074, -0.011045784173580318], gphase = [-0.8148890270369707, 2.910651779718185, 0.6061595527224816, -0.866896361773986, 2.978959739761571, 0.389449823323487, -0.9200529149001753, 3.042003639447468, 0.19086777256130605, -0.9730693796085267  …  -0.5845719433598235, -0.8020811950903951, 0.25019278693408326, -0.9442387897175247, 2.024847725056052, -0.5290350581915129, -0.9155932078752517, 0.30957465757960817, -1.1470666827470526, 2.0871896780327828]), (c = (params = [-0.5498584189890113 -0.19484289145814535 … -0.9077069338439421 0.05747367822846193; -1.354487379239275 -0.2091986952980046 … -0.3012650458531218 -1.5188481971774626; … ; 0.7734810216227409 0.6980885768675956 … -0.39763079146554703 0.18897020808620063; 1.2535985495609872 0.5745708256677979 … -0.9185895307323744 0.21080809906177214], hyperparams = 11.514372377669021), fg = 0.3239048517680847, σimg = 0.9802324204282875, lgamp = [-0.019863457976187986, 0.028179334624514624, -0.8601485485208558, -0.01252883535423731, -0.009762723982916771, 0.006876317277925958, -0.3255531196782044, 0.12462695299431686, 0.014681379345082065, -0.0114647145380579  …  -0.06774483157280944, -0.009301156007189145, -0.9433658654068834, 0.0038089761226709675, -0.021255580479178315, 0.01474442188134809, -0.06364734192634156, -0.0015639844974780666, -0.970371558323889, 0.00936702012064465], gphase = [-0.8126130480428638, 2.9374004624387706, 0.6160257569280557, -0.8673004417650652, 2.9928755975415857, 0.4058598661429304, -0.9291538707029982, 3.053804387709239, 0.22049883424141115, -0.9735498638444297  …  -0.5847908564661825, -0.780818014468177, 0.29259170131297935, -0.9125122182744378, 2.060307502818936, -0.5180569596984576, -0.9027896523187836, 0.3136206098458527, -1.1220430963557801, 2.099007265147497]), (c = (params = [0.20205176599962793 0.4711729282306738 … -0.5279231558161995 -1.258365917465042; -0.4899261753214675 -0.976277155193449 … -1.6817644486814876 -0.5348117114646103; … ; 0.621412503109798 0.7922297493094853 … -0.4096624761949101 -0.4673817076669111; -0.5787470473622292 0.5950141042522455 … -0.5528408742443063 -1.2232436863720857], hyperparams = 13.074582546523885), fg = 0.3495038839196328, σimg = 1.0057010662859314, lgamp = [-0.044158728230420066, 0.04814627203284011, -0.9122631801461057, 0.05086798422979616, -0.03779735597902874, 0.03825297819996588, -0.30498899472283, 0.15970477868430888, 0.008394253669737112, -0.009479999912346477  …  -0.07011550998504125, -0.0066408072686325184, -0.9446591973946715, -0.008927664051790742, -0.018174341443338, 0.017759831686544866, -0.07849163155751013, 0.006913270642491477, -0.9822968313882572, -0.00394421341580412], gphase = [-0.810798652730645, 2.9802770674531116, 0.6895267145537337, -0.8760645734125878, 3.04371114517513, 0.4805189825965584, -0.9273779617187733, 3.0992547502171317, 0.26993266278099604, -0.9814762063807003  …  -0.5834191515102093, -0.8048338329833898, 0.41011263001568077, -0.9440547146676492, 2.1942245053853218, -0.5178199470456849, -0.9513831514977593, 0.44834472846510737, -1.150168680380598, 2.233568713089146]), (c = (params = [0.4039786496213564 -0.3108566797262861 … -0.9477014558728241 -1.2015951666942473; -0.2962021989937289 -0.758304848553316 … -0.8548954519302099 -0.59895884014025; … ; -0.0666714778772098 1.35967768842461 … 0.461949847324379 -0.015160988706634912; -0.49952385907570795 1.1269045445968433 … -0.6980610573765552 -0.7026527508247369], hyperparams = 26.28636486808211), fg = 0.3602610208130456, σimg = 1.039457260147401, lgamp = [-0.03783018240756362, 0.03249421648485252, -0.8855431164315569, -0.014788082859595382, -0.029933491458037236, 0.034271547749173216, -0.3357082471855267, 0.11987313877134245, 0.01872796389031242, -0.020154403773246586  …  -0.05496073160775033, 0.0019430201616596563, -0.9452462538421372, 0.00018836136210293304, -0.01613015942263431, 0.022810331462362977, -0.07205916728125289, 0.006425624029712643, -0.9964097535915668, -0.0016096254181842348], gphase = [-0.8105200957024068, 3.005462195369246, 0.7033883918453059, -0.8765700248502516, 3.0572256791908865, 0.4832374035234097, -0.9279625217402634, 3.1185216406632845, 0.2818068878140039, -0.972462115815054  …  -0.5807623284130894, -0.7642114016946134, 0.38463768165147827, -0.9325552586670707, 2.1573149646804683, -0.5331897953656398, -0.8988772554187509, 0.4031313434655054, -1.1441919429148335, 2.185215318517793]), (c = (params = [-0.13986823921032793 -0.20638006554377405 … -0.673283452209818 -0.9189391992933048; -0.1679426774608821 -0.7403127664266815 … -1.1506086610961779 -0.38973123966147494; … ; -0.41289918013653204 1.4012614886785904 … 0.21253722270967523 -0.03294235499482614; -0.46546860198083656 0.8519846427367718 … -1.143633936947256 -0.5579874925318263], hyperparams = 31.61647541963119), fg = 0.3619820539833288, σimg = 1.0125739607255133, lgamp = [-0.03396784130642802, 0.03934779140445929, -0.9467792720732123, -0.014407953697521358, -0.025079758772899374, 0.0213975557275348, -0.3330666160654709, 0.12201894597587558, 0.020124795367675852, -0.01653453544520464  …  -0.07069077754048136, 0.0036920472747165828, -0.9223393446196562, 0.010867057405273758, -0.02169406402993132, 0.015548160009363788, -0.0820711667123497, 0.004483923974182538, -0.9779479575278254, -0.0038411278246805207], gphase = [-0.8158230587954017, 3.007714927802267, 0.7484873544314502, -0.8734826433544095, 3.0721941398758923, 0.5463664555234362, -0.9251381458897254, 3.1270830781395746, 0.32846867607601776, -0.9761491976986812  …  -0.581499452145762, -0.7551799003044143, 0.3845889605693079, -0.9488097466000395, 2.1481053791320472, -0.5141735689818789, -0.8942305790016778, 0.3916747059821577, -1.141501348273603, 2.1872504468667646])], NamedTuple{(:n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :max_hamiltonian_energy_error, :tree_depth, :numerical_error, :step_size, :nom_step_size, :is_adapt), Tuple{Int64, Bool, Float64, Float64, Float64, Float64, Float64, Int64, Bool, Float64, Float64, Bool}}[(n_steps = 1023, is_accept = 1, acceptance_rate = 0.9926747491491344, log_density = 2936.7531819127325, hamiltonian_energy = -2385.5792181617944, hamiltonian_energy_error = 0.007981392507190321, max_hamiltonian_energy_error = 0.010777337564377376, tree_depth = 10, numerical_error = 0, step_size = 0.0001, nom_step_size = 0.0001, is_adapt = 1), (n_steps = 2, is_accept = 1, acceptance_rate = 4.9807540489572726e-36, log_density = 2936.7531819127325, hamiltonian_energy = -2282.821522120973, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 4868.890043806264, tree_depth = 1, numerical_error = 1, step_size = 0.001554613274507503, nom_step_size = 0.001554613274507503, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.945324706341931, log_density = 2519.0275232065, hamiltonian_energy = -2242.0941316576927, hamiltonian_energy_error = 0.06968754294575774, max_hamiltonian_energy_error = 0.11197186170466011, tree_depth = 10, numerical_error = 0, step_size = 0.0003024688192109267, nom_step_size = 0.0003024688192109267, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9378287881594526, log_density = 2272.0634660145934, hamiltonian_energy = -1845.1244621462088, hamiltonian_energy_error = 0.0708804748987859, max_hamiltonian_energy_error = 0.19693302501241305, tree_depth = 10, numerical_error = 0, step_size = 0.0004354454109302323, nom_step_size = 0.0004354454109302323, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9426303142118211, log_density = 2043.6661441561291, hamiltonian_energy = -1637.6415542916973, hamiltonian_energy_error = 0.05292987483176148, max_hamiltonian_energy_error = 0.33484758618556043, tree_depth = 10, numerical_error = 0, step_size = 0.0007013302103589483, nom_step_size = 0.0007013302103589483, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.5831029077348284, log_density = 1941.1042444567508, hamiltonian_energy = -1387.596133133658, hamiltonian_energy_error = 0.27886217964578464, max_hamiltonian_energy_error = 5.315446556716552, tree_depth = 10, numerical_error = 0, step_size = 0.0012264256771442807, nom_step_size = 0.0012264256771442807, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.7173409536826532, log_density = 1730.9819388447277, hamiltonian_energy = -1214.102599076528, hamiltonian_energy_error = -0.12116575347681646, max_hamiltonian_energy_error = 1.1155713318116796, tree_depth = 10, numerical_error = 0, step_size = 0.0007397779351954428, nom_step_size = 0.0007397779351954428, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9888374172091614, log_density = 1582.1181688859601, hamiltonian_energy = -1054.9716473276458, hamiltonian_energy_error = -0.08955479414544243, max_hamiltonian_energy_error = -0.2131587863398181, tree_depth = 10, numerical_error = 0, step_size = 0.0006649385149941593, nom_step_size = 0.0006649385149941593, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.7507738439388639, log_density = 1486.9380566021875, hamiltonian_energy = -941.4099318259811, hamiltonian_energy_error = 0.1885727724602475, max_hamiltonian_energy_error = 1.2404816189209669, tree_depth = 10, numerical_error = 0, step_size = 0.0014029731171140093, nom_step_size = 0.0014029731171140093, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.3516238841336562, log_density = 1332.3154450997706, hamiltonian_energy = -806.8518080848426, hamiltonian_energy_error = -0.8569894776368301, max_hamiltonian_energy_error = 68.14355465847893, tree_depth = 10, numerical_error = 0, step_size = 0.0014087105469965796, nom_step_size = 0.0014087105469965796, is_adapt = 1)  …  (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9900525091224777, log_density = 1165.3243046444127, hamiltonian_energy = -483.4465032854986, hamiltonian_energy_error = -0.24403476602935825, max_hamiltonian_energy_error = -0.7098850386789763, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.5047681612837129, log_density = 1174.4805658508049, hamiltonian_energy = -497.5842676211719, hamiltonian_energy_error = 0.5121639649566987, max_hamiltonian_energy_error = 1.9068871304752975, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9189100852225911, log_density = 1170.8227603328987, hamiltonian_energy = -507.46876299272526, hamiltonian_energy_error = 0.13489553687190892, max_hamiltonian_energy_error = 0.5064829433727027, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.853695821339512, log_density = 1165.420657518465, hamiltonian_energy = -484.6850997323809, hamiltonian_energy_error = 0.10299811978995876, max_hamiltonian_energy_error = 0.716863308482516, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.7101194775599189, log_density = 1186.1628812926424, hamiltonian_energy = -514.5235664644492, hamiltonian_energy_error = 0.2568850060896466, max_hamiltonian_energy_error = 1.13879763477496, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.907756238424618, log_density = 1197.6618554406587, hamiltonian_energy = -550.4599001247861, hamiltonian_energy_error = -0.5413014000248495, max_hamiltonian_energy_error = -1.0007507955333494, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9527713830073807, log_density = 1185.344722772455, hamiltonian_energy = -499.06546783010174, hamiltonian_energy_error = -0.14236888138566428, max_hamiltonian_energy_error = -0.516972960276803, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.5669973938813859, log_density = 1175.3943334481555, hamiltonian_energy = -524.0784324268902, hamiltonian_energy_error = 0.6322921731963334, max_hamiltonian_energy_error = 1.543847192718772, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9191315822461662, log_density = 1155.9121566949916, hamiltonian_energy = -514.9068435606123, hamiltonian_energy_error = -0.3480050229477456, max_hamiltonian_energy_error = -1.6284785168920735, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9784005657940732, log_density = 1146.1817086598305, hamiltonian_energy = -477.52687708638985, hamiltonian_energy_error = -0.5014871772635843, max_hamiltonian_energy_error = -1.1656449489846636, tree_depth = 9, numerical_error = 0, step_size = 0.005928627584205218, nom_step_size = 0.005928627584205218, is_adapt = 0)])
Note

The above sampler will store the samples in memory, i.e. RAM. For large models this can lead to out-of-memory issues. To fix that you can include the keyword argument saveto = DiskStore() which periodically saves the samples to disk limiting memory useage. You can load the chain using load_table(diskout) where diskout is the object returned from sample. For more information please see ComradeAHMC.

Now we prune the adaptation phase

chain = chain[501:end]
Table with 5 columns and 200 rows:
      c                     fg        σimg      lgamp                 ⋯
    ┌──────────────────────────────────────────────────────────────────
 1  │ (params = [-0.26394…  0.265033  0.901203  [-0.0357336, 0.0400…  ⋯
 2  │ (params = [0.290632…  0.260263  0.892481  [-0.0498146, 0.0505…  ⋯
 3  │ (params = [-0.07242…  0.295311  1.00583   [-0.0162858, 0.0141…  ⋯
 4  │ (params = [0.756886…  0.33143   0.963831  [0.0113966, -0.0041…  ⋯
 5  │ (params = [-0.37625…  0.287716  0.968132  [-0.00773061, 0.006…  ⋯
 6  │ (params = [1.18051 …  0.288472  0.878922  [-0.0356138, 0.0410…  ⋯
 7  │ (params = [-0.38920…  0.282092  0.903824  [-0.0212832, 0.0235…  ⋯
 8  │ (params = [0.702054…  0.301829  0.904386  [-0.0217631, 0.0185…  ⋯
 9  │ (params = [0.301306…  0.282418  0.899126  [-0.0203189, 0.0180…  ⋯
 10 │ (params = [0.11909 …  0.31274   0.957545  [-0.0150051, 0.0193…  ⋯
 11 │ (params = [-0.18848…  0.319227  1.00951   [-0.0170908, 0.0127…  ⋯
 12 │ (params = [1.27048 …  0.307236  0.973562  [-0.00754161, 0.000…  ⋯
 13 │ (params = [0.287639…  0.313983  0.976292  [0.0319029, -0.0344…  ⋯
 14 │ (params = [2.08381 …  0.301793  1.01673   [-0.00243682, 0.002…  ⋯
 15 │ (params = [0.502946…  0.324628  0.931289  [-0.0157767, 0.0119…  ⋯
 16 │ (params = [0.659926…  0.322347  0.974949  [-0.00870677, 0.004…  ⋯
 17 │ (params = [1.32917 …  0.303714  0.975673  [-0.0198807, 0.0084…  ⋯
 ⋮  │          ⋮               ⋮         ⋮               ⋮            ⋱
Warning

This should be run for likely an order of magnitude more steps to properly estimate expectations of the posterior

Now that we have our posterior, we can put error bars on all of our plots above. Let's start by finding the mean and standard deviation of the gain phases

gphase  = hcat(chain.gphase...)
mgphase = mean(gphase, dims=2)
sgphase = std(gphase, dims=2)
104×1 Matrix{Float64}:
 0.003674223156644287
 2.6085714065905625
 0.19159212184291635
 0.0033325832945331375
 3.067095444346909
 0.19540423939707235
 0.0036014111214110506
 2.7315922290827785
 0.19873266353756117
 0.0031006754177965888
 ⋮
 0.10343630214875549
 0.1431796412186126
 0.08690490419438168
 0.14386570396982107
 0.004516595559860776
 0.10366323853588479
 0.1387975556099779
 0.08375840340159109
 0.1386788114152357

and now the gain amplitudes

gamp  = exp.(hcat(chain.lgamp...))
mgamp = mean(gamp, dims=2)
sgamp = std(gamp, dims=2)
129×1 Matrix{Float64}:
 0.017152158971260684
 0.018206962685153944
 0.02011435527796557
 0.035676289777485166
 0.015597552038124273
 0.015869256318966882
 0.027275527615870127
 0.037159329787201816
 0.01564563441260573
 0.015990870318815204
 ⋮
 0.00975430394443385
 0.01361866746494823
 0.009070312379444567
 0.01493523302857317
 0.015691659425243084
 0.02416431145859513
 0.009973455086300646
 0.013319277491710399
 0.009405742024434836

Now we can use the measurements package to automatically plot everything with error bars. First we create a caltable the same way but making sure all of our variables have errors attached to them.

using Measurements
gmeas_am = measurement.(mgamp, sgamp)
ctable_am = caltable(gcache, vec(gmeas_am)) # caltable expects gmeas_am to be a Vector
gmeas_ph = measurement.(mgphase, sgphase)
ctable_ph = caltable(gcachep, vec(gmeas_ph))
───────────┬────────────────────────────────────────────────────────────────────
      time │      AA            AP          AZ           JC            LM      ⋯
───────────┼────────────────────────────────────────────────────────────────────
 0.917±0.0 │ 1.0±0.0  -0.814±0.004     missing      missing       1.6±2.6    0 ⋯
 1.217±0.0 │ 1.0±0.0  -0.871±0.003     missing      missing    -0.006±3.1    0 ⋯
 1.517±0.0 │ 1.0±0.0  -0.929±0.004     missing      missing      -1.4±2.7    0 ⋯
 1.817±0.0 │ 1.0±0.0  -0.977±0.003     missing      missing      -2.0±2.3    0 ⋯
 2.117±0.0 │ 1.0±0.0  -1.002±0.004     missing      missing      -2.5±1.6    0 ⋯
  2.45±0.0 │ missing       1.0±0.0     missing      missing      2.6±0.12   0. ⋯
  2.75±0.0 │ 1.0±0.0  -1.046±0.007     missing      missing      -2.7±1.1   -0 ⋯
  3.05±0.0 │ 1.0±0.0  -1.133±0.003     missing      missing    -2.79±0.85   -0 ⋯
  3.35±0.0 │ 1.0±0.0  -1.152±0.003     missing      missing     -2.8±0.85   -0 ⋯
 3.683±0.0 │ 1.0±0.0  -1.163±0.005   2.46±0.11      missing      -2.6±1.4   -0 ⋯
 3.983±0.0 │ 1.0±0.0  -1.167±0.003   2.31±0.11      missing      -2.2±2.1  -1. ⋯
 4.283±0.0 │ 1.0±0.0   -1.16±0.003   2.14±0.11      missing      -1.1±2.9  -1. ⋯
 4.583±0.0 │ 1.0±0.0  -1.142±0.004   1.95±0.11    -1.1±0.13       1.7±2.6  -1. ⋯
 4.917±0.0 │ 1.0±0.0  -1.074±0.006   1.74±0.11   -0.88±0.14       2.6±1.3  -1. ⋯
 5.183±0.0 │ 1.0±0.0  -1.079±0.006   1.57±0.11   -0.72±0.14      2.76±0.1  -1. ⋯
  5.45±0.0 │ 1.0±0.0  -1.118±0.011   1.43±0.11   -0.56±0.14      2.53±0.1  -1. ⋯
     ⋮     │    ⋮          ⋮            ⋮            ⋮            ⋮            ⋱
───────────┴────────────────────────────────────────────────────────────────────
                                                    2 columns and 9 rows omitted

Now let's plot the phase curves

plot(ctable_ph, layout=(3,3), size=(600,500))
Example block output

and now the amplitude curves

plot(ctable_am, layout=(3,3), size=(600,500))
Example block output

Finally let's construct some representative image reconstructions.

samples = skymodel.(Ref(post), chain[begin:2:end])
imgs = intensitymap.(samples, fovx, fovy, 128,  128)

mimg = mean(imgs)
simg = std(imgs)
fig = CM.Figure(;resolution=(400, 400))
CM.image(fig[1,1], mimg,
                   axis=(xreversed=true, aspect=1, title="Mean Image"),
                   colormap=:afmhot)
CM.image(fig[1,2], simg./(max.(mimg, 1e-5)),
                   axis=(xreversed=true, aspect=1, title="1/SNR",),
                   colormap=:afmhot)
CM.image(fig[2,1], imgs[1],
                   axis=(xreversed=true, aspect=1,title="Draw 1"),
                   colormap=:afmhot)
CM.image(fig[2,2], imgs[end],
                   axis=(xreversed=true, aspect=1,title="Draw 2"),
                   colormap=:afmhot)
fig

Now let's check the residuals

p = plot();
for s in sample(chain, 10)
    residual!(p, vlbimodel(post, s), dvis)
end
p
Example block output

And viola, you have just finished making a preliminary image and instrument model reconstruction. In reality, you should run the sample step for many more MCMC steps to get a reliable estimate for the reconstructed image and instrument model parameters.

This page was generated using Literate.jl.